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I.  INTRODUCTION 


The  interpretation  of  thermal  infrared  (IR)  images  generated  from 
aircraft  or  satellite  sensor  data  requires  an  understanding  of  the  surface  and 
atmospheric  parameters  that  determine  the  thermal  response  of  the  ground  to 
the  incident  solar  insolation.  A  well-constituted  physical  model  that  accu¬ 
rately  simulates  the  energy  balance  at  the  earth-atmosphere  interface  is  an 
indispensable  analysis  tool  for  use  in  the  identification  of  surface  para¬ 


meters  that  cause  horizontal  signature  differences.  For  example,  if  indepen¬ 
dent  measurements  of  surface  shortwave  albedo  and  local  meteorological  vari¬ 
ables  are  available,  successive  iterations  of  a  thermal  modal  executed  over  a 
range  of  various  values  for  ground  conductivity  can  be  registered  with  a  back¬ 
ground  temperature  image  produced  from  IR  measurements  to  infer  the  horizontal 
distribution  of  the  thermal  conductivity  of  the  soil.  This  information  can  be 
used  to  determine  the  geologic  surface  materials  distributed  across  a  large 
image  scene.  In  a  similar  manner,  soil  moisture  and  anthropogenic  or  geo¬ 
thermal  heat  sources  can  be  inferred. 

For  use  in  background  modeling  applications,  a  modified  version  of  the 
boundary-layer  model  (HFLUX),  developed  by  Carlson  and  Boland*  for  simulation 
of  the  behavior  of  ground  temperature  over  complex  terrain,  is  being  imple¬ 
mented  at  The  Aerospace  Corporation.  The  model  can  be  used  in  various  field 
experiments  that  require  accurate  prediction  of  the  surface  radiometric 
temperature. 


II.  MODEL  STRUCTURE 


The  model  represents  the  physics  of  the  response  of  the  ground  to  the 
solar  forcing  by  the  use  of  Monin-Obukov  theory  with  an  implicit  K-type 
parameterization  for  the  eddy  fluxes  of  heat,  moisture,  and  momentum  in  the 
surface  layer.  The  primary  forcing  mechanism  is  the  net  radiation,  but  the 
partitioning  of  the  energy  balance  into  the  ground  heat  flux  and  the  latent 
and  sensible  heat  fluxes  into  the  atmosphere  is  determined  by  the  substrate 
and  atmospheric  variables.  Solutions  are  obtained  as  successive  equilibrium 
states  of  the  energy  balance  controlled  by  similarity  theory  in  the  surface 
layer  and  the  temperature  diffusion  equation  in  the  soil.  Typically,  the 
model  is  integrated  to  simulate  a  complete  diurnal  cycle  over  which  the  ter¬ 
rain  parameters  and  the  temperature  at  the  bottom  of  the  substrate  slab  are 
held  constant.  Horizontal  advection  is  not  included  in  the  model. 

A  schematic  view  of  HFLUX  is  presented  in  Fig.  1.  The  substrate  slab  is 
of  variable  depth  with  a  user-defined  number  of  levels,  permitting  coarse  or 
detailed  representations  of  the  subsurface  thermal  parameters.  The  surface 
air  layer  depth  is  fixed  at  50  m,  with  a  1-cm  transition  layer  at  the  earth's 
surface  that  contains  both  molecular  and  eddy  heat  conduction.*  Within  this 
transitional  sublayer,  the  fluxes  of  moisture  and  heat  are  assumed  to  be 
independent  of  the  stability.  A  well-mixed  layer  is  represented  in  bulk  form 
above  the  surface  layer  and  its  depth,  H,  is  governed  by  a  formulation  pre¬ 
sented  by  Tennekes^  in  which  H  grows  throughout  the  day  as  a  result  of  the 
surface  heat  flux  from  below  and  entrainment  of  air  with  a  higher  potential 
temperature  from  above.  The  atmosphere  above  the  surface  layer  is  modeled 
differently  at  night;  the  structure  of  the  nocturnal  boundary  layer  is  de¬ 
scribed  later  in  this  section. 

The  solution  sequence  begins  with  the  net  radiation  given  by  the  surface 
energy  balance 


The  exact  depth  of  the  transition  layer  is  Inconsequential  insofar  as  the 
surface  teaperature  determination  is  concerned. 
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(1) 
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RN  ”  Go  +  Ho  +  Eo 

where  GQ,  HQ,  and  EQ  represent  the  heat  flux  into  the  ground,  the  sensible 
heat  flux  into  the  atmosphere,  and  the  latent  heat  flux  into  the  atmosphere. 
All  quantities  are  positive  when  directed  away  from  the  surface.  The  net 
radiation  is  actually  calculated  from  the  radiative  fluxes  at  the  surface 


RN  “  S  +  Fd  "  Fu 


(2) 


where  S  is  the  solar  flux  at  the  surface  and  and  Fu  are  the  upward  and 
downward  terrestrial  fluxes,  respectively.  Figure  2  is  a  schematic  of  the 
energy  balance  that  shows  the  direction  of  the  fluxes  at  midday.  Total  down¬ 
dwelling  irradiance  absorbed  at  the  surface,  S,  is  calculated  from  a  one-layer 
radiative  transfer  model.  S  is  a  function  of  solar  geometry,  atmospheric 
transmission  coefficients,  and  ground  albedo,  AQ,  determined  as 


a  -  a  > 

S  =  S*  (3) 

o 

where  S*  and  X  are  general  transmission  relations  containing  the  solar  con¬ 
stant  and  coefficients  for  dust,  water  vapor,  air  molecules,  ozone,  and 
clouds.  S  can  be  calculated  for  sloping  terrain  of  arbitrary  orientation. 

The  upward  and  downward  directed  fluxes  of  longwave  terrestrial  radiation  are 
given  as 


F  •  VT/  <4) 

Fd  -  «a°T44  (5) 


where  a  is  the  Boltzmann  constant,  and  e  and  e,  are  the  emissivlties 

8  A 

ground  and  atmosphere.  Most  surface  materials  have  values  of  e  that 

g 

tween  0.9  and  1,  and  is  a  function  of  the  water  vapor  mixing  ratio 
surface  layer. 


for  the 
are  be- 
in  the 
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SW  *  +  IR  f  -  lfH  =  E  ♦  +  H  ♦  +  G  {  WHERE: 

SW  \  =  SOLAR  RADIATION 
IR  *  =  ATMOSPHERIC  INFRARED  RADIATION 
IR  f  =  SURFACE  INFRARED  RADIATION 
E  ♦  =  LATENT  HEAT  FLUX 

H  ♦  =  SENSIBLE  HEAT  FLUX 
G  I  =  GROUND  HEAT  FLUX 

EARTH-ATMOSPHERE 
INTERFACE 

Fig.  2.  Energy  Balance  at  Earth-Atmosphere  Interface 


a.-  ' 


.  r* 


The  surface  fluxes  of  sensible  and  latent  heat  are  parameterized  in  terms 
of  eddy  and  molecular  dif fusivities  as 


=  ~(cc 


+  pcpV  ft 


(6) 


E 


o 


+  pLFKn 

E  q 


\ 

;  3  z 


(7) 


where  Cg  and  Cw  are  the  molecular  diffusivities  for  heat  and  water  vapor,  q. 
Lg  is  the  latent  heat  of  vaporization,  and  Kg  and  are  the  eddy  diffusivi¬ 
ties  for  heat  and  water  vapor,  which  are  assumed  to  be  equal.  The  sensible 
and  evaporative  heat  fluxes  at  the  surface  are  calculated  by  integrating  Eqs. 
(6)  and  (7)  from  Z  =  0  to  the  top  of  the  surface  layer,  ZA 


PC. 

“o  ’  T*  (»„  *  V  W) 

Eo  ■  °-T-  M  -  ’„>  <9> 

where  qQS  is  the  surface  saturation  mixing  ratio,  qA  is  the  mixing  ratio  at 
ZA,  and  M  is  the  moisture  availability  parameter  described  by  Nappo.^  Values 
of  M  range  from  0  to  1  and  represent  a  fraction  of  the  potential  evaporation 
rate  for  a  saturated  surface.  As  such,  the  factor  M  accounts  for  the  reduc¬ 
tion  in  the  efficiency  of  evaporation  caused  by  subsaturation  of  the  surface. 
Evaporation  is  set  to  zero  at  night  when  qA  >  qos* 

In  Eq.  (8),  I  is  essentially  the  vertical  integral  of  the  inverse  of  the 
diffusivities 


1 


fZ  A  dz 

o'  kh  +  C./ocp 


1  1+  1  2 


(10) 
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This  integral  is  calculated  in  two  sections:  from  z  ■  0  to  ,  which  is.  the 
depth  of  the  transitional  layer  over  which  molecular  conduction  dominates  the 
transfer  process,  and  from  to  Zfi,  where  Cg  and  Cw  are  assumed  to  be  negli¬ 
gible  compared  with  pc  K  .  Thus,  turbulent  eddy  conduction  determines  the 

P  “ 

exchange  in  evaluating  the  expression  Ij.  The  term  ku*z/$^  is  substituted  for 
KH  where  u*  is  given  by 

“Ak 

u*  *  ~l - - -  (ID 

l  A  $  (z)/z(dz) 
l  tn 
o 

and  uA  is  the  wind  speed  at  ZA,  k  is  von  Kansan' s  constant  (0.4),  and  ZQ  is 

the  roughness  length.  The  parameter  u*  reflects  the  magnitude  of  mechanically 

produced  turbulence.  The  functional  forms  for  and  A  have  been  given  by 
,  H  M 

Panof  sky . 4 

The  ground  heat  flux  is  given  by  the  standard  conduction  equation 

A(T  -  T  )  (12) 

G  -  - 2- - — 

o  Az 


where  A  is  the  thermal  conductivity  of  the  substrate  and  T_j  is  the  tempera¬ 
ture  at  the  first  soil  level  below  the  surface.  The  transfer  of  heat  through 
the  soil  is  governed  by  the  diffusion  equation 


3T 

3t 


(13) 


where  <  is  the  thermal  diffusivity  of  the  substrate  and  is  equivalent  to 

A/C  ,  where  C  is  the  ground  heat  capacity. 

8  o 

The  diffusivity  and  conductivity  can  be  combined  to  form  a  parameter 
called  the  thermal  inertia,  P,  where 


P 


(14) 


12 


is  a  measure  of  the  rate  of  heat  transfer  at  the  interface  between  the  ground 
and  atmosphere.  Sensitivity  tests  indicate  that  whereas  model  solutions  de¬ 
pend  on  the  value  of  P  independently  of  the  choices  X  and  k,  X  and  ic  must  not 
be  chosen  completely  independent  of  one  another.  Values  of  X  and  ic  cited  in 
the  literature  appear  to  vary  systematically  for  a  wide  variety  of  real  mater¬ 
ials.  Accordingly,  20  pairs  of  X  and  k  based  on  values  reported  in  the  Manual 
of  Remote  Sensing^  and  by  Sellers**  have  been  assembled  to  produce  a  regression 
equation  whereby  specification  of  P  will  determine  uniquely  a  X,  k  pair. 

Thus,  it  was  found  that 


where 


X  =  -0.00013  +  0.0502P  +  1.21P2 


(15) 


< 


(16) 


Equation  (15)  explains  91%  of  the  variance  of  X  about  P  and  provides  an  em¬ 
pirical  result  that,  for  most  types  of  materials,  realistically  represents  the 
physical  relationship  of  X  and  tc  to  P.  The  regression  Eq.(15)  corresponds  to 
a  diverse  range  of  soil  types  (pumice  to  quartzite). 

Equations  (8),  (9),  and  (12)  can  be  combined  with  the  energy  balance  [Eq. 
(1)]  to  form  an  expression  for  the  sensible  heat  flux  at  the  surface 


Rjj  -  f(uA)M(qQS  -  qA>  -A 
Ho  '  1  +  DI 

in  which 

x<rzA  +  T0  -  T.,) 

A  *  - n - 


Azpc 

p 

where  T  is  the  dry  adiabatic  lapse  rate,  and  f(uA)  is  pLg/I. 


(17) 


(17a) 


(17b) 
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The  ground  temperature,  TQ,  is  determined  by  equating  the  right-hand 
sides  of  Eqs.  (1)  and  (2),  which  produces  the  quartic  equation 


where 


and 


AT  +  BT  +  C  -  0 
o  o 


A  ■  e  o 
g 


Az 


-1  4 

C  =  S  +  X— -  E  .oT .  -  H  -  E 

Az  A  A  o  o 


(18) 

(18a) 

(18b) 

(18c) 


At  each  time  step,  Newton's  technique  for  finding  real  zeros  of  a  polynomial 
is  iterated  until  Eq.  (18)  converges  to  the  correct  value  for  TQ.  At  night, 
in  the  nonturbulent  surface  layer,  H0,  EQ,  and  S  are  zero  so  that  Eq.  (18) 
expresses  a  balance  between  the  longwave  terrestrial  radiation  and  the  ground 
heat  flux.  As  such,  the  value  of  the  thermal  inertia,  P,  is  an  important 
determinant  in  the  nighttime  behavior  of  TQ. 

This  series  of  equations  is  cycled  through  for  convergence.  One  itera¬ 
tion  per  time  step  of  4  min  has  been  found  sufficient  to  ensure  accurate  solu¬ 
tions.  The  flow  chart.  Fig.  3,  schematically  traces  the  order  of  calculation; 
the  sequence  ends  with  the  determination  of  the  temperature  profile  in  the 
soil  at  the  next  time  step  from  which  the  model  integration  time  is  incremen¬ 
ted,  and  the  cycle  is  repeated.  The  temperature  profile  in  the  surface  layer 
and  substrate,  as  well  as  the  surface  radiative  and  turbulent  energy  fluxes, 
are  all  part  of  the  solution  set. 

The  nighttime  mode  of  calculation  requires  a  different  methodology,  be¬ 
cause  the  physics  driving  the  surface  layer  are  quite  dissimilar  from  those  of 
the  daytime.  During  the  day,  the  net  radiation  is  the  forcing  mechanism,  and 
the  heat  flux  component  is  a  function  of  the  moisture  availability  and  thermal 

14 
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inertia.  The  surface  layer  is  then  unstable  and  the  heat  flux  determines  the 
stability.  Near  sunset,  the  upward  heat  flux  vanishes  and  the  radiational 
cooling  causes  the  ground  temperature  to  decrease  rapidly;  subsequently,  the 
heat  flux  changes  sign  and  becomes  passively  dependent  on  the  lapse  rate  near 
the  ground.  To  account  for  this,  a  modified  form  of  a  scheme  proposed  by 
Blackadar^  has  been  incorporated  into  the  model.  The  maintenance  of  turbu¬ 
lence  at  night  is  calculated  as  a  function  of  the  bulk  Richardson  number,  B, 
in  the  surface  layer.  The  value  of  B  determines  the  form  of  the  profile 
equations  and  is  given  by 

Z  Z 

8  -  f  “T  [(6A  "  V  +  T*ln  Z~J  (19) 

W .  o 

A 

where  g  is  the  acceleration  caused  by  gravity,  6  is  an  average  temperature  in 
the  surface  layer,  0A  is  the  potential  temperature  at  ZA,  WA  is  the  magnitude 
of  the  wind  at  ZA,  and  0  ^  is  a  "shelter"  height  temperature  provided  by  a 
prognostic  equation  relating  radiational  convergence  and  turbulent  flux 
convergence 


-  0^  +  B 


H 

o 

PCPZ1 


(20) 


in  which  A  and  B  have  been  empirically  determined  as  8.3  *  10  ^s-*  and  0.2, 
respectively.  Monin-Obukov  scaling  is  used  to  obtain  T*,  u*,  and  the  surface 
heat  flux,  HQ.  Thus 


0A-01 


*  ln  VZA  -  *h 


(21) 


kW. 


*  ln  Z./Z  p 
1  o  Tm 


(22) 


H  -  -  kpc  u*T* 

n  n  *  * 


(23) 
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where  Z  is  the  roughness  lengthy  i|>u  and  \|>  are  the  nondimensional  profiles 
for  temperature  and  wind,  the  functional  forms  of  which  are  dependent  on 
stability,  and  is  the  magnitude  of  the  wind  vector  at  a  height  (50  m). 
The  vertical  profiles  for  temperature  and  wind  are  provided  through  the  inte¬ 
gration  of  the  U,  V  momentum  equations  and  the  thermodynamic  equation 


3U.  K  .  K  . 

TT  «  f(Vj  -  V  )  +  (04  +  j  -  U,)  -  (Uj  "  Di  -  i>  (24) 

Az  Az 


rr  -  -£<ui  -  V  +  -sLrJ-  (vi  +  i '  vi>  -  <5i  -  vi  - 1>  <25> 


^  +  I  -1,-^ 


where  f  is  the  Coriolis  parameter  and  Az  is  a  layer  depth  of  50  m.  Horizontal 
advection  has  not  been  Included  in  the  model.  The  eddy  exchange  coefficients 
for  momentum  and  heat,  KM  and  Kjj,  are  assumed  equal  in  the  stable  nocturnal 

Q 

boundary  layer  and  are  given  by  a  form  proposed  by  Mellor  and  Yamada  ,  with 
the  use  of  boundary  layer  data  and  second  order  closure 

Ki  “  1  Si  <27> 


where 


l(Ui  '  Ul-1)2  +  (Vi  "  Vi-l)2]1/2 


and  1  scales  as  kz  (taken  to  be  28  m)  in  the  surface  layer.  The  local 
Richardson  number  is  calculated  as 

.  .  1  (ei-9i-u 
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and  the  critical  Richardson  number  is  calculated  as  a  function  of  the  geo- 

Q 

strophic  wind  with  the  use  of  an  empirical  result  given  by  Vigeant 


0.5542  e-l°-2I29<V2  +  V2>1/21  +  0.2 


When  the  local  Richardson  number  exceeds  the  critical  value,  turbulent  ex¬ 
change  ceases  at  that  height  and  K  is  arbitrarily  set  equal  to  zero. 

The  lowest  layer  must  be  treated  differently,  because  there  is  no  turbu¬ 
lent  exchange  through  the  underlying  surface.  The  last  terms  in  Eqs.  (24) 

2  2 

through  (26)  must  be  replaced  by  -u^  U^/W^z,  -u^  N^/W^Az,  and  -H^/CpPAz, 
respectively. 

The  integration  is  performed  over  a  1000-m-deep  layer  consisting  of  20 
equally  spaced  intervals  (i  =  1,20)  50  m  apart  (Fig.  4)  and  is  begun  in  the 
late  afternoon  or  evening  when  the  surface  heat  flux  has  changed  sign.  For¬ 
ward  differencing  is  used,  and  a  time  step  of  120  sec  has  been  found  adequate 
to  ensure  computational  stability. 

Initially,  the  1000-m  layer  is  assumed  to  be  well  mixed,  and  the  vertical 
distribution  of  8  is  set  equal  to  a  constant  0A  at  the  time  of  the  stability 
reversal,  which  is  reasonable  for  late  afternoon  on  a  sunny  day.  The  vertical 
shear  of  the  geostrophic  wind  is  determined  from  an  analysis  of  local  surface 
fields  of  temperature  and  pressure;  usually  only  one  or  two  additional  wind 
observations  above  the  surface  layer  are  available.  To  fit  a  profile  for  the 
actual  winds  to  the  observations,  free  convection  scaling  is  assumed  to  govern 
the  wind  profile  before  the  upward  heat  flux  vanishes.  The  free  convection 
scaling  assumption  is  that  turbulence  in  the  lower  part  of  the  mixed  layer  is 
buoyantly  driven  (Ref.  10),  and  provides  a  relationship  for  the  wind  shear 
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which  results  in  a  velocity  defect  expression 


where  b  is  determined  empirically  and  z  is  measured  downward  from  z0bs  (where 
z  ■  0  and  U  ■  U^g)*  Equation  (32)  is  applied  to  00Z  wind  data  for  the  layer 
between  the  top  of  the  surface  layer  at  50  and  200  m;  above  200  m  the  wind  is 
assumed  initially  to  be  constant  up  to  1000  m.  The  complete  form  of  the  late 
afternoon  wind  and  temperature  profiles  is  represented  in  Fig.  5.  The  speci¬ 
fied  initial  wind  profile  may  possess  a  significant  ageostrophic  component, 
but  the  model  adjusts  smoothly  to  the  transition  during  the  first  1  or  2  hr  of 
integration  beyond  the  time  of  reversal  in  the  sign  of  the  heat  flux,  and  ul¬ 
timately  reaches  quasi-equilibrium  with  small  accelerations. 

The  value  of  the  bulk  Richardson  number  defines  three  stability  classi- 


fications  for  which  events 

in 

the  surface 

layer  are  different 

(I) 

B 

<  0 

Unstable 

(II) 

0 

<  B  <  0.2 

Stable 

(III) 

B 

>  0.2 

Nonturbulent 

In  cases  I  and  II,  the  surface  heat  flux  [given  by  Eq.  (23)]  is  positive  when 
the  surface  layer  is  unstable  (a  situation  that  rarely  occurs  at  night),  and 
negative  (downward  directed)  when  the  surface  layer  is  stable.  For  case  III, 
the  surface  layer  is  decoupled  from  events  above  50  m  because  of  the  strongly 


zero. 

A  typical  sequence  of  events  near  sunset  would  follow  a  progression 
through  the  three  stability  classifications  with  the  surface  layer,  initially 
unstable,  becoming  stable  (B  increasing).  As  the  ground  temperature  and  wind 
shear  stress  decrease  rapidly,  the  stratification  decouples  the  surface  layer 
from  the  atmosphere  above.  As  the  stress  decreases  and  the  turbulence  near 
the  ground  ceases  (B  >  0.2),  the  wind  above  accelerates,  increasing  the  verti¬ 
cal  shear.  If  this  shear  becomes  large  enough,  B  can  fall  below  0.2  and 
permit  a  turbulent  event  to  occur  with  a  downward  directed  exchange  of  heat, 
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which  will  act  to  elevate  the  ground  temperature  temporarily.  Such  turbulence 
may  recur  throughout  the  night,  Whereas  the  model  will  not  be  very  accurate 
in  timing  these  events,  numerical  tests  demonstrate  that  this  formulation  suc¬ 
cessfully  predicts  the  nights  that  are  likely  to  experience  turbulent  episodes. 
The  governing  parameter  for  these  events  is  the  geostrophic  wind  (assumed  con¬ 
stant  with  time).  On  nights  when  the  surface  geostrophic  wind  is  large,  tur¬ 
bulent  events  are  more  likely. 
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III.  MODEL  CAPABILITIES  DEVELOPED  AT  AEROSPACE 

A.  DUAL-MODE  CALCULATIONS 

In  many  remote-sensing  applications,  it  is  desirable  to  have  the  capacity 
to  simulate  the  thermal  behavior  of  a  relatively  small  target  area  surrounded 
by  terrain  with  very  different  surface  properties  (such  as  thermal  Inertia  and 
albedo),  e.g.,  an  asphalt  road  traversing  an  unvegetated  range  of  sandy  soil. 
If  the  target  surface  area  is  small  (less  than  a  few  hundred  square  meters), 
the  atmospheric  properties,  temperature  and  humidity,  overlying  the  target 
would  be  driven  by  the  atmospheric  fluxes  of  the  surrounding  terrain  rather 
than  the  temperature  and  humidity  of  the  small  target  (except  in  a  thin  layer 
above  the  target  where  the  eddy  scale  is  small  and  horizontal  advection  is 
minimal).  The  atmosphere  of  the  surrounding  terrain  thus  influences  the 
energy  balance  of  the  target  surface;  the  local  temperature  and  moisture 
anomalies  of  the  target  do  not  alter  the  large-scale  atmospheric  properties 
driven  by  the  energy  balance  of  the  surrounding  terrain. 

This  dual  mode  capacity  has  been  incorporated  into  a  version  of  the  model 
(HDUAL)  that,  in  one  recent  application,  was  effective  in  predicting  the  ther¬ 
mal  signature  of  a  concrete  pad  surrounded  on  all  sides  by  a  composite  sand- 
granite  soil.  The  target  and  terrain  are  modeled  to  possess  separate  (and 
often  quite  dissimilar)  values  for  parameters  such  as  thermal  inertia,  albedo, 
emissivity,  and  moisture  availability. 

B.  SUBSURFACE  MODEL 

In  its  original  form  HFLUX  used  an  explicit  differencing  scheme  to  model 
the  diffusion,  Eq.  (13),  that  governs  the  transfer  of  heat  through  the  sub¬ 
strate.  The  stability  criterion  for  this  differencing  scheme  limits  the  ver¬ 
tical  resolution  permissible  in  the  substrate,  and  the  thermal  properties  were 
restricted  to  be  constant  with  depth.  In  many  applications  this  is  adequate, 
but  in  the  current  version  of  HFLUX,  the  former  differencing  method  has  been 
replaced  by  the  implicit  Crank-Nicholson  scheme  to  permit  better  resolution  in 
the  substrate  and  allow  for  depth-dependent  values  of  thermal  inertia  and 
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multiple  layers  (of  different  materials).  With  this  procedure,  it  is  now 
possible  to  model  a  complex  substrate  matrix,  and  applications  that  require 
fine  resolution  near  the  surface  are  easily  handled.  This  submodel  also 
permits  the  inclusion  of  multiple  boundary  conditions  that  are  necessary  in 
special  simulations. 

C.  DYNAMIC  UPDATING 

HFLUX  now  has  the  capacity  to  accept  dynamic  updating.  Meteorological 
field  data  measured  at  5-min  intervals  have  been  reduced  and  placed  on  a  data 
file.  These  data  are  read  sequentially  into  the  model  at  each  time  step  and 
replace  the  customary  theoretical  calculations  of  solar  insolation,  atmo¬ 
spheric  radiation,  and  subsurface  temperature  within  HFLUX  with  actual  weather 
station  measurements.  With  dynamic  updating,  model  calculations  are  nudged 
with  measured  field  data  toward  the  actual  thermal  state  of  the  ground  and 
atmosphere,  and  this  results  in  a  more  accurate  prediction  of  ground  tempera¬ 
ture  and  the  surface  energy  balance.  Also,  the  handling  of  atmospheric  tran¬ 
sients  such  as  cloudiness  is  far  superior  to  that  of  the  empirical  relations 
that  would  otherwise  be  used  to  simulate  their  presence.  Tests  comparing 
updated  and  nonupdated  simulations  on  days  with  significant  cloudiness  and 
variable  weather  have  conclusively  demonstrated  the  utility  of  dynamic 
updating. 
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IV.  FUTURE  DEVELOPMENTS 


A.  SOIL  MOISTURE 

Sensitivity  tests  indicate  that  the  amplitude  of  the  ground  temperature 
wave  is  most  responsive  to  changes  in  the  surface  moisture  availability  and 
thermal  inertia.  The  value  of  these  parameters  is  closely  related  to  the 
distribution  of  moisture  in  the  upper  layer  of  the  soil  (~  10  cm),  and  the 
capacity  to  model  water  transport  in  the  soil  would  be  of  great  value.  The 
rate  of  liquid  water  and  water  vapor  movement  in  a  soil  is  a  strong  function 
of  the  soil  composition  (e.g.,  sand  and  clay),  and  the  determination  of  the 
coefficients  in  the  transport  equations  that  are  unique  to  each  soil  type  will 
require  ouch  more  extensive  field  work.  This  constitutes  a  fundamental 
problem  in  the  creation  of  a  soil  moisture  submodel.  The  coefficients  exist 
in  the  literature  for  a  limited  number  of  soils  that  have  received 
attention. ^ 

Deardorff  has  proposed  a  predictive  equation  for  the  surface  moisture 
derived  from  the  conservation  equation  for  moisture 


5Wg 
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(Wg  -  w2) 


in  which  Wg  is  the  surface  moisture  content,  W2  is  the  vertically  averaged 
value  of  moisture  over  the  upper  25  cm  of  soil,  and  Tj  is  the  diurnal  period. 

and  C2  are  coefficients  that  vary  with  moisture  content  and  soil  type.  If 
Ci  and  C2  were  known  for  various  soils,  the  variation  of  moisture  availability 
with  time  as  a  function  of  evaporation  and  water  transport  could  be  calculated 
for  a  broad  class  of  soil  types.  When  the  moisture  content  is  known  for  the 
top  10  cm,  the  soil  thermal  properties  can  be  determined  as  a  function  of 
moisture  content.  In  a  saturated  soil,  water  fills  the  gaps  between  soil 
granules,  increasing  the  bulk  thermal  inertia,  whereas  air  occupies  these  gaps 
in  a  dry  soil,  which  results  in  relatively  small  values  for  P.  Pratt  and 
Ellyett^  have  discussed  in  detail  the  determination  of  thermal  inertia  as  a 
function  of  moisture  content. 
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A  version  of  Deardorff's  submodel  soon  will  be  included  in  HFLUX  to  simu¬ 
late  the  behavior  of  the  surface  soil  moisture  and  its  effect  on  the  thermal 
properties  of  the  soil. 

B.  ENERGY  BUDGET  OVER  WATER 

Another  extension  of  HFLUX  that  may  prove  useful  in  future  applications 
will  include  the  capacity  to  determine  the  surface  energy  balance  over  water. 

A  current  version  (HPUDL)  calculates  the  surface  temperature  and  evaporation 
time  scale  for  a  small  puddle  of  water  over  an  asphalt  surface. 

C.  MODEL  CONSOLIDATION 

A  longer  terra  objective  is  to  construct  a  very  generalized  form  of  the 
model  that  would  permit  the  simulation  of  thermal  signatures  for  a  complex  set 
of  terrain  features  such  as  a  matrix  of  buildings,  water  bodies,  vegetation, 
and  homogeneous  or  natural  materials,  such  as  concrete  and  sand.  Successive 
application  of  the  model  for  each  of  the  above  scenes  would  produce  a  thermal 
image  that  would  be  very  useful  in  scene  identification  and  image  decorrela¬ 
tion  techniques. 
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